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Abstract 

We present a novel sparse modeling approach to non-rigid shape match- 
ing using only the ability to detect repeatable regions. As the input to 
our algorithm, we are given only two sets of regions in two shapes; no 
descriptors are provided so the correspondence between the regions is 
not know, nor we know how many regions correspond in the two shapes. 
We show that even with such scarce information, it is possible to estab- 
lish very accurate correspondence between the shapes by using methods 
from the field of sparse modeling, being this, the first non-trivial use of 
sparse models in shape correspondence. We formulate the problem of per- 
muted sparse coding, in which we solve simultaneously for an unknown 
permutation ordering the regions on two shapes and for an unknown cor- 
respondence in functional representation. We also propose a robust vari- 
ant capable of handling incomplete matches. Numerically, the problem is 
solved efficiently by alternating the solution of a linear assignment and a 
sparse coding problem. The proposed methods are evaluated qualitatively 
and quantitatively on standard benchmarks containing both synthetic and 
scanned objects. 



1 Introduction 



Matching of deformable shapes is a notoriously difficult problem playing an im- 
portant role in many applications KZHCO10 . Unlike rigid matching where the 



correspondence can be parametrized by a small number of parameters (rotation 
and translation of one shape w.r.t. the other CM91 , BM92]), non-rigid match- 
ing typically uses point-wise representation of correspondence, which results in 
the number of degrees of freedom growing exponentially with the number of 
matched points. 

Non-rigid correspondence methods try to find correspondence by minimizing 
some structure distortion. The structures can be point-wise (local descriptors 



ZBVH09||SOG09||GB AL09| ASC 1 II ) , p air- wise (distances |ETC||MS05l|BBK06 
BBK^ToJ), or higher order [ZWW+ 10 



In order to make the matching problem computationally feasible, it is crucial 
to reduce the size of the search space 



TBW+11 



Most methods use a combi- 
nation of point- and pair-wise structure matching in order to achieve this, and 
typically consist of three main components: feature detection, feature descrip- 
tion, and regularization. Given two shapes, a feature detector allows to find a 
set of landmarks (points or regions) that are repeatable, i.e., appear (possibly 
with some inaccuracy) on both shapes. Feature descriptor then assigns to each 
feature a vector capturing some local geometric properties of the shape; very 
often, the two processes are combined into a single one. Using the descriptors, 
landmarks on two shapes can be matched (it has been shown OMMG10 that 
under some conditions, correct landmark matching fully determines the intrin- 
sic correspondence between the shapes). Such a matching reduces the search 
space size to points with similar descriptors. However, since the matching uses 
only local information, such correspondence can be noisy, and some kind of 
regularization based on higher-order information is needed to rule out bad or 
inconsistent correspondences. This information is also used to establish the 
correspondence between the rest of the points on the shapes. Often, the pro- 
cess is applied hierarchically, restricting the candidate matches to points in the 



proximity of the landmarks SY12 



Computer graphics and geometry processing literature contains a plethora of 
approaches for each of the aforementioned components. Feature detection meth- 
ods try to locate stable points or regions [DMAMS10 LBB11 that are invariant 
under isometric deformations and robust to noise. Popular feature descrip- 
tors include the heat kernel signature (HKS) |SOG09 |GBAL09 , wave kernel 
signature (WKS) |ASC11 , global point signature (GPS) [Rus07 or methods 



adopted from the domain of image analysis [ZBVH09| . As regularization, pair- 



wise structures such as geodesic MS05 , BBK06 or diffusion distances BBK+ 10 
and higher-order structures ZWW + 10 have been used. 

Alternatively, there have been several attempts to represent correspondences 
with a small set of parameters. Elad and Kimmel EK01 used multidimensional 
scaling (MDS)-type methods to embed the intrinsic structure of the shapes into a 
low-dimensional Euclidean space, posing the problem of non-rigid matching as a 
problem of rigid matching of the corresponding embeddings ("canonical forms"). 
Mateus et al. MHK+08 used spectral embeddings instead of MDS. Lipman and 
Funkhouser LF09 embedded the shapes into a disk by means of conformal maps 
and represented the correspondence as a Mobius transformation. 
More recently, Ovsjanikov et al. 



OBCS+12 



introduced the functional rep- 
resentation of correspondences, allowing to perform a "calculus" of correspon- 
dences. In this approach, correspondence is modeled as a correspondence be- 
tween functions on two shapes rather than points, and can be compactly repre- 
sented in the Laplace-Beltrami eigenbasis 



Lev06 as a matrix of coefficients of 



decomposition of the basis functions of the first shape in the basis of the second 
one. In this paper, we will be relying upon this latter representation. 
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1.1 Main contribution 



The main practical contribution of this paper is an approach for finding dense 
intrinsic correspondence between near-isometric shapes with very little known 
information: we only assume to be able to detect regions in two shapes in a 
repeatable enough way (i.e., that at least some regions in one shape correspond 
accurately enough to some other regions in another shape). No region descrip- 
tors are given, so the correspondence of the regions is unknown. The assumption 



of near-isometry assures that in the functional representation of OBCS + 12 , the 



unknown correspondence can be represented as a sparse matrix. The assump- 
tion of repeatable regions implies that there exists some unknown permutation 
that orders the regions according to their correspondence. 

We formulate the problem of permuted sparse coding, in which we simulta- 
neously look for the permutation and the correspondence, thereby introducing 
the very successful area of sparse modeling into efficient and state-of-the-art 
shape correspondence. We note that with the permutation fixed, our problem 
becomes the standard sparse coding problem; having the correspondence fixed, 
the problem becomes a linear assignment. This allows efficient numerical so- 
lution by alternating the two aforementioned problems and employing efficient 
solvers that exist for both. 

Our method relies on a pretty common assumption that the shapes are 
nearly-isometric (though our experimental results show our approach still works 
even when departing from this assumption) , and out of all methods we are aware 
of, it uses perhaps the scarcest amount of data to establish dense correspondence 
between the shapes. For example, sandard region detectors with high repeata- 
bility such as LBB11 are sufficient. 

Compared to recent techniques for region-wise shape matching (see, e.g. 



GF09 VKTS+11 HKG11 PBB11 ), our approach has several important prac- 



tical advantages: First, we do not use any feature descriptor. Second, most 
region-wise correspondence approaches require an additional step of extending 
the correspondence between matched regions to the rest of the points. 

The rest of the paper is organized as follows. In Section 2, we overview the 
functional representation of correspondences, allowing to work with correspon- 
dences as algebraic structures, and state the main notions in sparse modeling. 
In Section 3, we formulate our problem of permuted sparse coding for estab- 
lishing correspondence from a set of repeatable regions given in unknown order. 
We then extend the problem to the general setting where the region detection 
process is not perfectly repeatable. In Section 4, we describe the numerical 
optimization used to solve our permuted sparse coding problem. Experimental 
results are shown in Section 5. Finally, Section 6 discusses the limitations and 
possible extensions of the proposed framework and concludes the paper. 
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Figure 1: Near isometric shape correspondence as a sparse modeling problem 
(see details in text): Indicator functions of repeatable regions on two shapes are 
detected and represented as matrices of coefficients A and B in the corresponding 
orthonormal harmonic bases $ and 4'. When the regions are brought into 
correspondence, the point-to-point correspondence between the shapes can be 
encoded by an approximately diagonal matrix C. In the proposed procedure 
termed as permuted sparse coding, we solve IIB = AC + O simultaneously for 
an approximately diagonal C and the permutation II bringing the indicator 
functions into correspondence. To cope with imperfectly matching regions, we 
relax the surjectivity of the permutation and absorb the mismatches into a row- 
wise sparse outlier matrix O. For visualization purposes, the coloring of the 
regions is consistent as after the application of the permutation. Correspondence 
is shown between a synthetic TOSCA and scanned SCAPE shape. 




Figure 2: Top row: representation of different maps between the two shapes 
using the matrix C. Shown left-to-right are the ideal correspondence, a symmet- 
ric correspondence, and a random correspondence. Bottom row: representation 
of the same correspondences as point-to-point maps. Note that the farther is 
the correspondence from an isometry, the less diagonal is the matrix C. 



2 Background 

2.1 Functional representation of correspondences 

The direct representation of correspondences as maps between two non- Euclidean 
spaces limits the range of tools that can be employed for correspondence com- 
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Figure 3: Dense point-to-point correspondences obtained between the left 
TOSCA human shape and its approximate isometries. Corresponding points 
are marked with consistent colors. The average correspondence distortion is 
depicted in units of the reference shape diameter. The highest distortions are 
obtained on the non-isometric joints, but do not exceed 6% of the diameter. 



putation due to the lack of an algebraic structure. In this paper, we rely on the 
functional representation of correspondences introduced in OBCS + 12 , which 



overcomes this limitation. In what follows, we briefly review the main idea of 
such functional representations. 

Let X and Y be two shapes, modeled as compact smooth Riemannian man- 
ifolds, related by a bijective correspondence t : X — > Y. Then, for any real 
function / : X —> R, we can construct a corresponding function g : Y — >• M 
as g = f ot . The correspondence t uniquely defines a mapping between two 
function spaces T : J-(X, M) — > .F(Y,R), where J-(X,M.) denotes the space of 
real functions on X. Such a representation is linear, since for every pair of 
functions /i,/2 and scalars ai,a 2 , 



T(ai/i + a 2 / 2 ) = (ai/i + a 2 h) ° t 1 

= ai/i o t^ 1 + a 2 / 2 o t^ 1 
= a 1 T(f 1 ) + a 2 T(f 2 ). 

Assuming that X is equipped with a basis {<fii}i>i, any / : X 
represented as 

/ = aj^i 
i>i 



(1) 



can be 



(2) 



with the di being some representation coefficients (in case of an orthonormal 
basis, a, = (f,<f>i); in the general case, the coefficients are found by projecting 
the function / on the bi-orthonormal basis). Due to the linearity of T, 



T (/) = T £ = E a * T (^) ( 3 ) 

\i>l J i>l 

If the shape Y is further equipped with a basis {ipj}j>i, then for every i there 
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exists coefficients c^- such that 

T(&)=5>^;, (4) 

and we can write 

T(f) = ^2 a i c ij^j- ( 5 ) 

i,J>l 

Let us now assume finite sampling of X and Y, with m samples (for sim- 
plicity, we assume that the shapes are sampled at the same number of samples 
m. The extension to the case with a different number of samples is straightfor- 
ward). The bases are represented as the m x n matrices <& and $ containing, 
respectively, n discretized functions fa and ipj as the columns. The functions 
/ and g = T(f) can now be represented as n-dimensional vectors f = $a and 
g = \l>b with the coefficients a and b. Using this notation, Equation ([5| can be 
rewritten as 

*b = T(*r) = *C T a; (6) 
since 9 is invertible, this simply means that 

b T = a T C. (7) 

Thus, the nx n matrix C fully encodes the linear map T between the functional 
spaces, and contains the coordinates in the basis ^! of the mapped elements of 
the basis <&. 



2.2 Point-to-point correspondence 

Point-to-point correspondences assume that each point i on I corresponds to 
some point j on Y. In functional representation, this is equiv alent to ha ving 
C that makes each row of \&C T coincide with some row of $ OBCS+12 . In 
applications requiring point-to-point correspondence, given some C, it can be 
converted into a point-to-point correspondence by thinking of $ and ^ as n- 
dimensional points clouds, and orthogonal matrix C as a rigid alignment trans- 
formation between them. This procedure is equivalent to iterative closest point 



(ICP) in n dimensions OBCS + 12 , initialized with the given Co: first, for each 
row i of \&Co T , find the closest row j* in $ (this operation can be performed effi- 
ciently using approximate nearest neighbor algorithms) . Then, find orthonormal 
C minimizing J^i ~ *C T |j2 and set Co = C. This operation is repeated 
until convergence and can be performed efficiently over all the vertexes of X 
and Y using approximate nearest neighbor algorithms. 

2.3 Sparse modeling 

One of the main tools that will be used in this paper are sparse models. In 
what follows, we give a very brief overview of this vast field, and refer the reader 



to ElalO for a comprehensive treatise. The central assertion of sparse modeling 



6 



is that many families of signals (and later operations as here introduced) can be 
represented as a sparse linear combination in an appropriate domain, usually 
referred to as the dictionary. This can be written as x « Dz, where x denotes the 
signal, D the dictionary, and z the sparse vector of representation coefficients. 
The dictionary is often selected to be overcomplete, i.e., with more columns than 
rows. 

Finding the representation of a signal x in a given dictionary D is usually 
referred to as sparse representation pursuit or sparse coding. Among the variety 



of pursuit methods, we will focus on the so-called Lasso formulation Tib96 
that finds z as the solution to the unconstrained convex program 

mm ||x - Dz||| + A||z||i. (8) 

z 

The hrst term is the data fitting term, while the second term involving l\ norm, 
||z||i = \z\ \ + . . . + \ z n \, promotes a sparse solution; the parameter A controls the 
relative importance of the latter. Proximal splitting methods Nes07 are among 
the most efficient and most frequently used numerical tools to solve problem ([8]); 
in Section |4j we present a variant of the proximal splitting algorithms for the 
solution of the pursuit problem arising in shape correspondence as detailed in 
the sequel. 

In some cases, signals not admitting the simplistic model of element-wise 
sparsity still manifest more intricate types of structured sparsity. In structured 
sparse models, the non-zero elements of z come in groups or, more generally, in 
hierarchies of groups. A common class of structured pursuit problems can be 
formulated as convex programs of the form 

min||x-Dz||^ + A||z||i, 2 , (9) 

z 

where the £12 norm, ||z|| li2 = IN1II2 + ' + ll z /c||2: assumes that the vector z 
is decomposed into k non-overlapping sub- vectors Zj, and promotes group- wise 
sparse solutions (i.e., the solution will have a small number of groups with non- 
zero coefficients, but the sub- vectors representing each such non-zero group will 
be dense). 

While structured sparse models enforce group structure of each representa- 
tion vector independently, it is often useful to consider the structure shared by 
multiple vectors. Collaborative sparse models operate on data matrices X, in 
which each column corresponds to a data vector, and assert that the patterns of 
non-zero coefficients are shared across the corresponding representation vectors, 
Z. This is achieved by solving a pursuit problem of the form 

min||X-DZ||| + A||Z|| 2)1) (10) 
z 

where the first term involving the Frobenius norm serves as the data fitting 
term, and the second term with the £2,1 norm promotes row-wise sparsity of 

the solution. The £2,1 norm is defined as 1 1 Z 1 1 2, a = \\ Z I\\2 H + INmlbj where 

zj denotes the i-th row of Z (note the difference from the £1,2 column-wise 
counterpart!). 
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Figure 4: Dense point-to-point correspondences obtained between the left 
SCAPE human shape and various other poses. Corresponding points are marked 
with consistent colors. 




Figure 5: First row: point-to-point correspondences obtained between different 
non-isometric shapes: male and female (left); two strongly non-isometric defor- 
mations of the dog shape from the TOSCA set (middle); TOSCA and SCAPE 
human shapes (right). Second row: Point-to-point correspondences obtained 
between SHREC shapes undergoing nearly isometric deformations and (from 
left to right) spike noise, Gaussian noise, and topological noise in the form of 
large and small holes. 



In this paper, we use formulate the shape correspondence problem using a 
sparse model, and use sparse modeling tools to efficiently solve it. 

3 Sparse modeling of correspondences 

In case the shapes X and Y are isometric and the corresponding Laplace- 
Bcltrami operators have simple spectra (no eigenvalues with multiplicity greater 
than one), the harmonic bases (Laplacian eigenfunctions) have a compatible be- 
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havior, ^ = T((f>i) such that = ±5y. Choosing the discretized eigenfunctions 
of the Laplace-Beltrami operator as $ and \t causes every low-distortion cor- 
respondence being represented by a nearly diagonal, and therefore very sparse, 
matrix C. 

In practice, due to lack of perfect isometry and numerical inaccuracies, the 
diagonal structure of C is manifested for the first eigenfunctions corresponding to 
the small eigenvalues (low frequencies) , and is gradually lost with the increase of 
the frequency (Figure [2]). However, a correspondence with low metric distortion 
will usually be represented by a sparse C. We use this property to formulate 
the computation of correspondences in terms of a sparse representation pursuit 
problem. 

Let us assume to have some region (or feature) detection process that given a 
shape X produces a collection of functions {fi : X — > M.} based on the intrinsic 
properties of the shape only. For example, the /j's can be indicator functions 



of the maximally stable components (regions) of the shape LBB11 . Since the 
process is intrinsic, given a nearly isometric deformation Y or X, it should 
produce a collection of similar functions {gj : Y — > H.}. 

To simplify the presentation, let us assume that the process is perfectly 
repeatable in the sense that it finds q functions on X and Y, such that for every 
fi there exists a g j = fi°t related by the unknown correspondence t. We stress 
that the ordering of the /j's and g^s is unknown, i.e., we do not know to which 
gj in Y a fi in X correspond. This ordering can be expressed by an unknown 
q x q permutation matrix II (in Section 3.2, we consider the more general case 
when the number of functions detected on X and Y can be different, i.e., II is 
non-square) . 

Representing the functions in the bases on each shape, we have fi = and 
gj = *&bj. Since each pair of corresponding fi and gj shall satisfy Q, we can 
write 

IIB = AC, (11) 

where A and B are the q x n matrices containing, respectively, &J and bj as 
the rows, and 7T,j = 1 if &i corresponds to hj and zero otherwise. 

3.1 Permuted sparse coding 



Note that in relation (11 1, both n and C arc unknown, and solving for them is 
a highly ill-posed problem. However, by recalling that the correspondence we 
are looking for should be represented by a nearly-diagonal C, we formulate the 
following problem 

minillllB-ACIII + AIIWoCIl!, (12) 

where the minimum is sought over n x n matrices C (capturing the correspon- 
dence t between the shapes in the functional representation) and q x q permu- 
tations II (capturing the correspondence between the detected regions on the 
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shapes). The first term containing the Frobenius norm can be interpreted as the 
data term, while the second term, containing the weighted l\ norm promotes a 
sparse C; denotes element-wise multiplication, and the non- negative param- 
eter A determines the relative importance of the penalty. Small weights Wij in 
W are assigned close to the diagonal, while larger weights are selected for the 
off-diagonal elements. This promotes diagonal solutions. 



The solution of ( 12 1 can be obtained using alternating minimization iterating 



over C with fixed II, and II with fixed C. Note that with fixed II, we can denote 



B' = HB and reduce problem ( 12 1 to 



min^lB'-AClH + A||W©C||i, (13) 

which resembles the Lasso problem frequently employed for the pursuit of sparse 
representations. On the other hand, when C is fixed, we set A' = AC, reducing 
the optimization objective to 

||IIB-A'||| = (14) 
tr (B T n T nB) - 2tr (B T II T A') + tr (A' T A') . 

Since II is a permutation matrix, II T II = I, and the only non-constant term 



remaining in the objective is the second linear term. Problem ( 12 ) thus becomes 



maxtr (II T E) , (15) 

where E = A'B T = ACB T and the maximization is performed over permutation 
matrices. To make it practically solvable, we allow II to be a double-stochastic 
matrix, which yields the following linear assignment problem: 

max vec(E) T vec(II) s.t. < J^i. \ (16) 
n>o II 1 = 1. 



We refer to problem (12 1 as to permuted sparse coding, and propose to solve 



it by alternating the solution of the standard sparse coding problem ( 13 1 and 



the solution of the linear assignment problem (16 1. The sparsity constraint has 
a regularization effect on this, otherwise extremely ill-posed, problem, and the 
following strong property holds: 



Proposition 1. The process alternating subproblems and {16) converges 
to a global minimizer of the permuted sparse coding problem Xli 



Due to lack of space, we provide the proof in the Appendix. This result 
means, among other, that despite the relaxation of the permutation matrix to 



a double-stochastic matrix in the assignment subproblem (16), we are actually 



recovering a true permutation matrix. This follows from the total unimodularity 



of the constraints in ( 16 ) 
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3.2 Robust permuted sparse coding 

So far, we have assumed the existence of a bijective, albeit unknown, correspon- 
dence between the fa's and the g/s. In practice, the process detecting these 
functions (e.g., stable regions) is often not perfectly repeatable. In what fol- 
lows, we will make a more realistic assumption that q functions fi are detected 
on X, and r functions gj detected on Y (without loss of generality, q < r), such 
that some /j's have no counterpart gj, and vice versa. This partial correspon- 
dence can be described by a q x r partial permutation matrix II in which now 
some columns and rows may vanish. 

Let us assume that s < q fi's have corresponding <7j's. This means that 
there is no correspondence between r — s rows of B and q — s rows of A, and 
the relation IIB sa AC holds only for an unknown subset of its rows. The 
mismatched rows of B can be ignored by letting some columns of II vanish, 
meaning that the correspondence is no more surjective. This can be achieved 



by relaxing the equality constraint II 1 = 1 in (16) replacing it with II 1 < 1. 



However, dropping injectivity as well and relaxing III = 1 to III < 1 would 
result in the trivial solution II = 0. To overcome this difficulty, we demand 
every row of A to have a matching row in B, and absorb the r — s mismatches 



in a row-sparse q x n outlier matrix O that we add to the data term of ( 12 ) 
This results in the following problem 

1 

mm - 
C,Oli 2 



Dim -nIIB - AG - 0||§. + A||W © C||i + /*l|0|| 2 ,i> ( 17 ) 



which we refer to as robust permuted sparse coding. The last term involves the 
£2,1 norm 

r 

\\Oh,i = £>?|| 3 , (18) 

i=l 

which can be thought of as the £\ norm of the vector of the £2 norms of the 
rows of O. The £2,1 norm promotes row-wise sparsity, allowing to absorb the 
errors in the data term corresponding to the rows of A having no corresponding 
rows in B; the parameter \i > controls the amount of regularization. The qx r 
matrix II is searched over all injective correspondences. 



As before, problem (17) is split into two sub-problems, one with the fixed 
permutation n, 

min^llB' -AC-0||p + A||W 0C|| 1 + H|0|| 2 ,i, (19) 

with B' = IIB, and the other one with the fixed C, 

max vec(E) T vec(n) s.t. \ ^T n = 3 n (20) 
n>o ( 11 1 < 1, 

with E = (AC)B T . Note that an injective correspondence is relaxed into a row- 
wise stochastic and column-wise sub-stochastic matrix II. Proposition 1 simply 
extends to the robust formulation as well. 
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4 Numerical solution 



The solution of the robust permuted sparse coding problem (17) is reduced 
to alternating two relatively standard optimization problems, and there exist 
many readily available efficient numerical tools to solve them. For the sake of 
completeness, we provide a concise description of the involved numerics. 



Problem ( 20 1 , being a simple linear assignment problem, is solved using the 
Hungarian algorithm. As an alternative, linear programming can be employed. 
To reduce the search space size, we disallow certain impossible permutations 
such as those relating regions with very distinct sizes. 



In order to solve (191, we use the family of forward-backward splitting al- 
gorithms Nes07 designed for solving unconstrained optimization problems in 
which the cost function can be split into the sum of two terms, 

min/i^x) +h 2 (x), (21) 

one, h±, convex and diffcrcntiable with an a-Lipschitz continuous gradient and 
another, hi, convex extended real valued and possibly non-smooth. Clearly, 



problem ( 19 1 falls in this category. 

The forward-backward splitting method with fixed constant step defines a 
series of iterates, {x 



X fe+1 =P„fc„ X fe - -^h„(„k 



l Vhi{^)\, (22) 



a 

where 

P Q/l2 (x) =argmin ||u - x||| + ah 2 (vL) (23) 

denotes the proximal operator of h 2 . Many alternatives have been studied in 
the literature to improve the convergence rate of forward-backward splitting 



algorithms BT09 Nes07 . Accelerated versions reach quadratic convergence 
rates (the best possible for the class of first order methods). The discussion of 
theses methods is beyond of the scope of this paper. 

In our case, the objective comprises a quadratic function hi = ||B'— AC— 0||| 
and the non-smooth function h 2 = A||W©C||i+/i||0||2,i. The proximal operator 
splits into two operators, one in C and another one in O, both having closed 
forms. The proximal operator corresponding to the ii norm term is given by 
the weighted soft threshold function 

Pi(C) = max||C|-^w|©sign(C), (24) 

where the absolute value and the sign functions are applied element-wise. The 
i-th row of the proximal operator corresponding to the i 2 i norm term is given 
by 



(P 2 (0)) 4 = max{||o^|| 2 -^}^. (25) 
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input : Data B',A; parameters A, /i; step size a. 

output: Sparse matrix O and row-wise sparse outlier matrix O 

Initialize 0° = B' and C° = 0. 

for k=l,2,. . . , until convergence do 

C fe+1 = Pi ((I - ^A T A)C fe - iA T (O fc - B')) 
O fc+1 = P 2 ((1 - %)O k ~ i(AC fc - B')) 
end 



Algorithm 1: Forward-backward splitting method for the solution of (19 1 



The gradient of the quadratic data term with respect to C and O is given 
straightforwardly by 

Vc/ii = A T AC + A t O A T B' 

V /ii = O + AC-B'. (26) 

The Lipschitz constant of the gradient determining the step size is bounded by 
the maximum eigenvalue 

/ a t A A T \ 

a < A max ( j \ . (27) 

Plugging the above expressions together into p2| yields the forward-backward 
splitting optimization summarized in Algorithm 111 

5 Experimental results 

In order to evaluate our approach, we performed several experiments on data 



from the TOSCA BBK08], SHRECT1 EB and SCAPE ASK+05 datasets 



The TOSCA set contains high-quality (10K-50K vertices) synthetic triangular 
meshes of humans and animals in different poses with known ground truth cor- 
respondences between them. SHREC'll contains meshes from the TOSCA set 
undergoing simulated transformations. The SCAPE set contains high-resolution 
(12K vertices) scans of a real human in different poses. 

For each pair of shapes we calculated the MSERs using 6-10 eigenfunctions 
and selected regions with areas of at least 5-10% of the total shape area, resulting 
in about 5 — 15 detected regions (see Figure [T]). These parameters were selected 
empirically for our data sets. 

The segments of each shape were projected onto 20 eigenfunctions and the 
corresponding C matrix was calculated by solving the sparse coding subprob- 
lem ( 19 ) using an accelerated variant of the method described in Section |4j 



The linear assignment subproblem (16) was solved using the Hungarian method 



Kuh55 . We initialized the permutation matrix with n = ^H T , and the cor- 
respondence matrix with C = 0. We observed a rapid convergence of the al- 
ternating minimization procedure in one or two iterations (see Figure [6] where 
for visualization purposed, II was initialized to identity). We found that the 
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Figure 6: Outer iterations of robust permuted sparse coding alternating the 



solution of the sparse representation purusit problem ( 19 1 with the linear as 



signment problem (20). Three iterations, shown left-to-right, are required to 



achieve convergence. Depicted are the permutation matrix II (first row), the 
correspondence matrix C (second row), and the outlier matrix O (last row). 
The resulting point-to-point correspondence and the correspondence matrix C 
refined using the ICP as described in Section |2.2| are shown in the rightmost 
column. 



method consistently converged to the same solution regardless of the initializa- 
tion. Finally, after convergence of the alternating minimization, the resulting C 



was refined using the method described in Section 2.2 



The robustness of the method is demonstrated in Figures [3]-[5j correct cor- 
respondences are computed even when the shapes undergo non-isometric defor- 
mations and are contaminated by geometric or topological noise. Figure [7] shows 
a quantitative evaluation and comparison of our algorithm to other correspon- 
dence algorithms on the SCAPE data set. The evaluation was performed using 
the code and data from 



KLF11 



Our method outperforms existing methods 
while making less assumption and working only with intrinsic information. 

The code used in the experiments was implemented in Matalb with parts 
written in C. The approximate nearest neighbor search in the ICP refinement 
step was accelerated using the FLANN library. The experiments were run on 
a 2.4GHz Intel Xeon CPU. End-to-end execution time varied from 10 to 50 
seconds, with the detailed breakdown summarized in Table [T] 



6 Discussion and Conclusion 

In this paper, we posed the problem of finding intrinsic correspondence between 
near-isometric deformable shapes as a problem of sparse modeling. Given only 
two set of regions in the two shapes with unknown correspondence, we are 
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Figure 7: Quantitative evaluation of the proposed shape correspondence algo- 
rithm and its comparison to other correspondence algorithms on the SCAPE 



shapes using the evaluation protocol from KLF11 



able to infer a dense correspondence between the shapes from two assumptions: 
that at least some of the regions in the two sets are corresponding; and that 
the shapes are nearly-isometric. The latter assumption implies that in func- 
tional representation in harmonic bases the unknown correspondence between 
the shapes is modeled as a sparse nearly-diagonal matrix; the former assumption 
implies that there exists an unknown permutation that reorders the regions in 
corresponding order. To find both the permutation and the correspondence, we 
formulate the novel permuted sparse coding problem and propose its efficient 
solution. An additional sparse coding term addressing outliers is added to the 
model for handling partial matching and formulated as the robust permuted 
sparse coding. 

To the best of our knowledge, among other dense correspondence techniques, 
our method relics on the smallest amount of information (the ability to find some 
repeatable regions) and quite generic assumption (near-isometric shapes). In 
particular, it allows us to use only a region detector without a feature descriptor 
to find a high-quality correspondence between two shapes. 



We note that, as in |OBCS + 12 , we explicitly assume that the shapes are 
nearly isometric, and that their Laplacians have simple spectrum. This assump- 
tion assures that the Laplacian eigenbases <& and $ have a compatible behavior, 
and as a result C has a nearly-diagonal structure. If we try to relax the restric- 
tion on multiplicity, C will still be sparse, but with unknown sparse structure. 
We can can still use our problem in this setting, imposing a different sparsity 
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Table 1: Average runtime (in seconds) as a function of the shape size for different 
stages in the proposed method: Basis - harmonic basis computation; MSER - region 
detection; Opt. - alternating minimization procedure; Ref. - ICP-based refinement 
and point-to-point correspondence computation; Tot. - total runtime. 



Vertices 


Basis 


MSER 


Opt. 


Ref. 


Tot. 


5K 


0.53 


0.61 


7.80 


1.41 


10.35 


10K 


0.99 


1.32 


7.91 


2.70 


12.92 


20K 


2.03 


3.58 


7.91 


5.52 


19.04 


50K 


5.57 


14.23 


7.85 


13.99 


41.64 



constraint on C. 

Relaxing the assumptions even more, we can depart from the near-isometric 
model, e.g. considering applications where one wishes to match shapes with 
roughly similar geometry but very different details (such as a horse and an 
elephant). In such a generic setting, the Laplacian eigenbases may differ dra- 
matically, and thus C have a non-sparse structure. It is possible to incorporate 
the bases <& and *tf as variables into our problem, and in addition to finding the 
permutation II and correspondence C find also the bases in which C will have 
a diagonal structure. This problem is akin to dictionary learning used in the 
sparse modeling literature. In future research, we will study such a generaliza- 
tion of our framework in hope to find correspondences between non-isometric 
shapes. Another possible generalization of our problem is for finding correspon- 



dence between collections of shapes NBCW + 11 



KLF11 



Finally, it worthwhile noting that the novel structured sparse modeling 
techniques introduced in |SBS12| provide an alternative to the optimization- 
based pursuit by replacing the iterative proximal algorithm with a learned 
fixed-complexity feed-forward network. Approaching shape correspondence as 
a learning problem from this perspective seems a very attractive future research 
direction. 



Appendix A — Proof of Proposition 1 

The permuted sparse coding problem 

mm JllIIB-ACIII+AllWoCII, (28) 

where P(g) denotes the space of q x q permutation matrices, is non-convex since 
the feasible set is non-convex. However, by relaxing P(q) to the bigger space of 
qx q double-stochastic matrices, S(q) = {IL € R q + Xq : III = n T l = 1} D P(g), 
we obtain the problem 



1 

Cei9 x ' ,fies(g) 2 



dud -||nB-AC||| + A||W0C|| (29) 
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that is easily shown to be convex. Fixing one of the variables at a time, the 
problem can be split into two subproblems: the sparse coding problem 

mini||B'-AC||f + A||W0C||i (30) 

C 2^ 

with B' = IIB, and the linear assignment problem 

max vec(E) T vec(n) (31) 
nes((j) 

with E = A'B T = ACB T . These two subproblems can be viewed as minimizing 



the objective function of (29) with respect to two blocks of coordinates, C and 



II. The minimization process alternating between the solution of (30) and (31) 
can be therefore regarded as block-coordinate descent. 



We use an instance of Theorem 4.1 in TseOl stating that block-coordinate 
descent is guaranteed to converge to a global minimizer of a non-differcntiablc 
convex function. The conditions of the theorem are satisfied by the objective 



and the constraints of ( 29 1 . Note that since we do not prove strict convexity of 



the latter objective, a global minimizer is not necessarily unique. 

Let us now have a closer look at the linear assignment problem (31) that 
can be cast as the linear program 

min e T 7r s.t. Qir = 1, (32) 
weu.f 

where the variable vector ir = vec(II) is the column stack of the assignment 
matrix, the cost vector is given by e = — vec(E), and the constraint matrix can 
be expressed using the Kronecker product notation as 

q = ( ) . <33) 

Here, I stands for the q x q identity matrix, and 1 T for the 1 x q vector of ones. 
In Lemma 1 below, we prove that Q is totally unimodular, i.e., any of its square 
submatrices has the determinant of or ±1. This property guarantees that the 
linear program has a solution with integer coordinates. Since E>(q)nZ qxq = ¥(q), 



this guarantees that the linear assignment problem (31) is actually equivalent 
to the binary assignment problem 

max vec(E) T vec(n). (34) 
neP(q) 

Combined this result with the convergence of the block-coordinate descent to 
a global minimizer of ( |29[ ), we can guarantee that in the obtained minimizer 
II actually belongs to P(q). This implies that the block-coordinate descent 



converges to a global minimizer of ( 28 1 



Lemma 1. The matrix Q in is totally unimodular. 



17 



Proof. The matrix Q can be constructed as the sub-matrix of a bigger matrix 



Q = ( V ) ® ( ]T ) ■ ( 35 ) 

Hence, total unimodularity of Q implies total unimodularity of Q. Since total 
unimodularity is preserved by the Kronecker product, it is sufficient to show 
that each of the two Kronecker factors are totally unimodular. We will limit the 
discussion to the second factor; very similar arguments apply to the first one. 
The matrix 

R = ( / T ) (36) 

comprises two components: the identity matrix I and the row vector 1 T , both 
of which are totally unimodular. For each k x k square submatrix R' of R, 
we distinguish between the following three cases: R' containing only elements 
of I; R' containing only elements of 1 T ; and R' containing elements of both 
components. In the two former cases, detR' € {0, ±1} since I and 1 T are 
totally unimodular. In the latter case, the submatrix has the form 



R' = O O , (37) 




where the to x to (to < k) identity matrix I is surrounded by zeros and concate- 
nated to a row of ones on the bottom. For m < k — 1, the submatrix contains 
at least one row of zeros and therefore detR' = 0. For m = k — 1, R' assumes 
the form 

R' = (J /t ) , (38) 

where I is the (k— 1) x (k — 1) identity matrix, is the (k — 1) x 1 vector of zeros, 
and 1 T is the 1 x (k — 1) vector of ones. Using the properties of the determinant, 
we obtain 

detR' = dct(01 T -I) = (39) 
Hence, R is totally unimodular. □ 
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